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Abstract 



We have performed deformed Hartree-Fock+BCS calculations with the Skyrme SIII force 
for the ground states of even-even nuclei with 2 < Z < 114 and N ranging from outside the 
proton drip line to beyond the experimental frontier in the neutron-rich side. We obtained 
^ ' spatially localized solutions for 1029 nuclei, together with the second minima for 758 nuclei. 

. The single-particle wavefunctions are expressed in a three-dimensional Cartesian-mesh rep- 

resentation, which is suitable to describe nucleon skins, halos, and exotic shapes as well as 
properties of ordinary stable nuclei. After explaining some of the practical procedures of the 
' calculations, we compare the resulting nuclear masses with experimental data and the pre- 

\^ . dictions of other models. We also discuss the quadrupole (m=0, 2) and hexadecapole {m—0, 

0^ ' 2, 4) deformations, the skin thicknesses, the halo radii, and the energy difference between 

the oblate and the prolate solutions. Our results can be obtained via computer network. 

^ ■ 1 Introduction 

It is one of the most important goals of the nuclear theory to reproduce and predict the nuclear 
ground-state energy and other properties globally in the nuclear chart within a single framework. 
^ ' A variety of theoretical models have been introduced for this goal. The most elaborate works 

have been done in the framework of the finite-range droplet model with a microscopic shell cor- 
rection (FRDM), whose latest result was given by Moller et al. Another extensive calculation 
was carried out by Aboussir et al. Q in the extended Thomas-Fermi plus Strutinsky integral 
method (ETFSI). The former as well as the latter methods can be regarded as approximations 
to the Hartree-Fock (HF) equation. The straight-forward solutions of the equation including 
deformation require long computation time for global calculations even with present computers. 
Such global results are not yet available to the public as far as we know. In this paper, we 
report on the results of our extensive HF-I-BCS calculation with the Skyrme SIII force for 1029 
even-even nuclei having atomic number ranging from 2 to 114 and neutron number from outside 
the proton drip line to beyond the experimental neutron-rich frontier. 

We use a HF-I-BCS code EV8 in which the single-particle wavefunctions are expressed in a 
three-dimensional Cartesian-mesh representation, while most of the other methods for deformed 
nuclei express the single-particle wavefunctions by the expansion in a harmonic oscillator basis. 
Let us mention three advantages of the mesh representation. First, it is capable of treating 
nucleon skins and halos. In contrast, they cannot be described efficiently in the oscillator- 
basis expansion because the asymptotic form of wavefunctions far from the nuclear surface is 
determined by the basis. Second, one can treat exotic (e.g., high-multipole) shapes and large 
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(e.g., super and hyper) deformations without preparing a basis specific to each shape. Third, 
the saturation property of density makes the atomic nucleus a very suitable object to the mesh 



representation, as we discuss in section 2.1 



This paper is mainly concerned with the proton-rich nuclei than the neutron-rich ones. We 
thoroughly explore the proton-rich side, even to beyond the proton drip line by several neutrons, 
while restricting the investigations in the neutron-rich side within the experimental frontier plus 
a few neutrons. This is because the HF+BCS method used in this paper cannot correctly 
include the coupling in the pairing channel with the continuum, which can be influential in 
nuclei near the neutron drip line[^, |5|. The HF-t-BCS method including the coupling with the 
continuum gives rise to unphysical neutron gas surrounding the nucleus. The solution of this 
problem requires the Hartree-Fock-Bogolyubov method, with which ordinary and pair densities 
are spatially localized when the Fermi level is at a negative energy Q. On the other hand, for 
nuclei near the proton drip line, the HF+BCS method is still applicable in practice to obtain 
localized solutions because the Coulomb barrier confines the wavefunctions of low-and-positive 
as well as negative energy levels. 

The contents of this paper are as follows. In the next section, we survey the features of 
the three-dimensional Cartesian-mesh representation, discuss the choice of the interaction, and 
explain the points newly developed in the present paper for extensive calculations. Among 
these points are the determination of the pairing force strengths and the acceleration of the 
convergence to the solution. 

In section ^, we compare the atomic masses obtained from our calculations with those in 
the latest experimental mass table[^]. For this purpose, we consider a correction for the error 
in the total binding energy due to the finite mesh size: A greatly high precision is necessary in 
calculating the binding energy compared with other quantities like moments. 

Section |^ treats the deformation. We compare the electric axial quadrupole moments with 
those deduced from the experimental B(E2)| dataQ. We also discuss the deformation parame- 
ters 020, 0^22) 0^40) 0,42, and a44, which we define for the HF+BCS solutions in terms of multipole 
moments, and compare them with those of the FRDM. The difference of shapes between pro- 
tons and neutrons as well as the difference of the energies between the prolate and the oblate 
solutions are also discussed. We choose a shape coexisting nucleus 4oZr4o to illustrate the strong 
dependence on the force parameters of the landscape of the potential energy curve along the 
axial-quadrupole-deformation path. 

The proton and the neutron skins are the subject of section ^. We discuss the thickness, the 
anisotropy, and the relation with the halo radius. 

In the last section, the contents of this paper is summarized and an instruction to obtain 
our results electronically via computer network is given. 



2 The set up of the calculations 
2.1 The basis 

The main feature of the HF+BCS code EV8 used in this paper is the three-dimensional Cartesian- 
mesh representation: Each single-particle wavefunction ip{x, y, z) is defined in a rectangular box 
{—\Lx < X < \Lx, —\Ly < y < \Lij-, —\Lz < z < \Lz) with its values il^ijk at cubic mesh 
points, {xi,yj, Zk) = (i — |, j — ^, A; — ^)a, where i, j, and k take on integers. In this study, the 
mesh size a is set to 1 fm, while the size of the box is = Ly = 26 fm, = 28 fm for Z < 82 
and Lx = Ly = 28 fm, = 30 fm for Z > 82. The nucleus is placed at the center of the box. 

We impose a symmetry with respect to reflections in x-y, y-z, and z-x planes (the point group 
D2h)- This symmetry allows triaxial solutions, although all of our solutions have eventually 
turned out axial and stable against 7-deformation. On the other hand, the symmetry prohibits 
odd-multipole deformations, which may not be negligible in some actinide nuclei. According to 
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the calculations with the FRDM[^, the nucleus ^ggRai34 has the largest octupole deformation 
(/?3 = 0.15, the energy gain due to the octupole deformation is —1.4 MeV), while except in 
the neighborhood of this nucleus the octupole deformation occurs only in odd-A and odd-odd 
nuclei. Incidentally, both in light |^] and heavy nuclei]^], the octupole deformation is likely to be 
enhanced by the procedure of the variation after parity projection. 

One might wonder that a mesh size of 1 fm were too large to describe the abrupt change 
of density at nuclear surface. It was demonstrated in Ref. |^, however, that a mesh size a=l 
fm can produce enough accurate results for several spherical nuclei with mass below ^'^^Pb. 
We did a similar test of accuracy for a deformed actinide nucleus ^94Pui46 and found that the 
relative errors of the quadrupole moment and the total energy are 0.4% and 0.5%, respectively. 
(The method of extrapolation to a —> is explained in section |3.1| .) This order of accuracy 
is higher than necessary for the quadrupole moment, while it is not for the energy to make 
comparison with experiments. Considering that the root-mean-square (r.m.s.) deviation of the 
atomic masses of recent mass formulae is ~ 0.5 MeV, the desirable precision is of the order of 0.1 
MeV, which is only 0.005% of the total binding energy of ^^''Pu. Therefore, the binding energy, 
but not the other quantities, has to be corrected for the effect of the finite mesh size, which is 



done in section 3.1 



The origin of this unexpectedly high accuracy with apparently coarse meshes has been ex- 
plained by Baye and IIeenen [p^ . The equation to determine {ipijk} is usually derived through a 
discrete approximation to the Schrodinger equation. They presented an alternative point of view, 
in which they introduced a set of orthogonal basis functions fijk{x,y, z) such that {ipijk} are 
the coefficients to expand ip{x, y, z) in this basis. (In this point of view, the equation for {ipijk} 
is determined uniquely from the variational principle.) This basis can be unitary-transformed 
to plane-wave basis with Ife^l < vr/a (k = x, y, z), which suggests that the atomic nucleus is a 
very suitable physical object to apply the mesh representation because the saturation property 
of nuclear matter guarantees the suppression of large-momentum components in wavefunctions 
from the view point of the Thomas-Fermi approximation. 

To enjoy this high accuracy, the method to determine {V'ijfc} must be in accordance with the 
view point of Baye and Heenen. Exact variational treatment with the plane-wave basis requires, 
however, long computation time and diminish the simplicity of the Cartesian-mesh representa- 
tion. The code EV8 is designed to emulate the plane- wave expansion method, though it is based 
on the discrete approximations, by choosing the appropriate orders of approximation formulae 
for derivatives (the 7- and 9-point formulae for the first and second derivatives, respectively) 
and integrals (the mid-point formula). 

It is worth mentioning that a new formulation of the mesh representation in terms of collo- 



cation basis splines is being developed recently by Chinn et al.|ll] 



2.2 The interaction 

For the HF (mean-field) part of the interaction, the code uses the Skyrme interaction |jl^, p!^ , 
which is a zero-range force with the lowest order momentum dependences to emulate the finite- 
range effects, a density dependence to reproduce the saturation of nuclear matter density, and 
a spin-orbit coupling term. 



The relation between the Skyrme-HF model and the relativistic mean-field model|14| has 
been discussed by many authors (see, e.g., Ref. |]l5|), which has invoked arguments on the 
density dependence (the ratio of the isoscalar to the isovector density dependences) of the spin- 



orbit term|16|, 17, 19|. However, as these arguments are not yet conclusive, we rather like to 
use the old but well-examined standard form of the Skyrme force in this paper. 

Among the many parameter sets proposed for the Skyrme force|Q, ^, 2C, 21, 2^, ^3, 24, 25, 



we choose the SIII||2^. Its validity has been examined in many nuclear structure calculations. 
In particular, it produces single-particle spectra in good agreement with experiment. It also 
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reproduces fairly well the N — Z dependence of the binding energy p7|] compared with other 



widely-used parameter sets of SGII|^ and SkM*[^. Although its incompressibility is said to be 
too large, it is not a serious drawback because elaborate fittings of the parameters of the Skyrme 
force and the FRDM[Q] have shown that different assumptions for the incompressibility lead 
to practically the same quality of fittings to nuclear masses. 

The force SkSC4[^ was determined through the most extensive fitting to nuclear mass data 
among the Skyrme forces. However, the fitting was done in the ETFSI scheme, which produces 
unnegligibly different energy from that of the HF method. Because of this disadvantage, we do 
not choose the force SkSC4 in this paper. 

2.3 The pairing 

For the interaction in the pairing channel, we employ a seniority force V^^jj., whose pair-scattering 
matrix elements are defined as a constant multiplied by cutoff factors depending on the single- 
particle energy e^, 

H^p^irliJ) = -Grfr{e^)fr{eJ), (1) 

where r signifies proton or neutron. We assume the following form for the cutoff function, 



/.(6) = a + exp-^-^V e{el-e), (2) 



-1/2 

0.5 MeV / 

which contains two cutoff energies, and ej: The former is for a smooth cutoff necessary to 
stabilize the iterative procedures to solve the HF equation, while the latter is for a sharp cutoff 
preventing occupation of spatially unlocalized single-particle states. The cutoff energies take on 
the following values, 

el = AJjF + 5 MeV, = + 2.3 MeV, (3) 

where A|jp is the Fermi level defined as the average of the highest occupied level and the lowest 
vacant level in the HF (or normal) state. 

In the BCS treatment of nuclei far from the /^-stability line, one has to take care so that the 
continuum states are not occupied, which give rise to unphysical nucleon gas extending over the 
entire box. In our calculations, for neutrons, if the right-hand side of the second of Eqs. (^) is 
positive, e° is replaced by zero (for Z < 82) or ej? is replaced by —2.3 MeV (for Z > 82). For 
protons, instead of lowering the sharp cutoff energy e^, we modify the proton potential outside 
the Coulomb barrier (in the imaginary-time evolution operator, not in the evaluation of the 
energy) so that the potential is higher than e^. This prevents the tunneling through the barrier 
and makes the proton single-particle wavefunctions whose energy is within the pairing active 
interval (e < ej?) spatially localized. 

In order to treat a wide range of nuclei on a single footing, we need a prescription to determine 
the strength Gr for each nucleus. For this purpose, we have developed a method based on the 
continuous spectrum approximation. We solve the following particle-number and gap equations, 

Nr = I {1- ^ ^ }D{e)de, (4) 

/(e-A,)2 + /,(e)2A2^ 




Ar = ^ / mde, (5) 

e-A;)2 + /,(6)2A2 

which use the semiclassical single-particle level density D{e) (obtained in the Thomas-Fermi 
approximation) instead of the discrete spectrum of the single-particle HF hamiltonian. One 
notes that l)(e) does not have the fluctuation due to the shell effects. 
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Eq. determines the Fermi level A,-, while Eq. (|5|) tells the force strength G-^ when it is 
combined with the empirical formula for the pairing gap = 12^-1/2 MeV. For hght nuclei, 
the resulting Gr becomes apparently too strong. We replace Gr with 0.6 MeV when it exceeds 
0.6 MeV. The results are discussed in section 3.3. 



2.4 The method of solution 

Because the single-particle basis of the mesh representation is huge, it takes unmanageably long 
computation time to solve the HF+BCS equation with the usual method of iterative diago- 
nalization of the single-particle HE hamiltonian /irf- Instead, the code EV8 employs a much 



more efficient method of imaginary-time evolution 1 28], in which the evolution operator for a 
small time interval At is repeatedly operated on each single-particle wavefunction to decrease 
its expectation value of the energy. After each evolution, the single-particle wavefunctions are 
orthogonalized in the Gram-Schmidt method from low to high energy states. As the initial 
wavefunctions, we utilize either the eigenstates of the Nilsson model or the solutions for the 
neighboring nuclei. 

The imaginary-time evolution operator exp(— /i~^/iHF^i) is expanded to the first order in 
At as 1 — fi~^h-ii-p^t- To this order, the imaginary-time evolution method is equivalent to the 



gradient iteration method] 29 1. Consequently, to obtain the minimum energy state rather than 
the maximum one, it must hold that 1 — /i~^emaxAt > —1, where emax is the maximum single- 
particle eigenenergy, which is approximately the free kinetic energy for k = S^/^vra^^. With a=l 
fm, it follows At < 2.1 x 10^^^ sec (The actual upper bound for At is larger than this estimation 
because the discrete approximation underestimates the kinetic energy for very high momentum 
states). We use At = 1.5 x 10~^^ sec in this paper. 

We regard that the wavefunction is converged to a HE+BCS solution when the following 
four criteria are satisfied. 

i) The energy spreading of single-particle states are smaller than 0.1 MeV. Practically, we 
require that the third largest value of 

(Aei)' = ((^l/iHpK) - (i|/iHFK)') • min (2vl l) (6) 

should be less than (0.1 MeV)^, where \i) represents the ith single-particle state and vf its BCS 
occupation probability. 

ii) During the last 75 steps of evolution, the variation in the total energy is < 1.5 keV for 
Z < 82 and < 5 keV for Z > 82. 

iii) During the same period, the variation in (x^)t-, {y^)T, and {z^)t are < 0.5 fm^ for Z < 82 
and < 0.005 x ^r^NrA?^^ for Z > 82, where r stands for all the protons or all the neutrons 
and ro = 1.2 fm. 

iv) The axial quadrupole deformation parameter 6 is close to the convergent value. We 



adopt a definition 6 = 3 {Qz) / 4(r )|3C], where Qz = 2z — x — y . The convergent value 



(5pred5 together with its uncertainty, Adp^ed, are predicted using the values of 5 during the last 
200 stepsQ. The condition is expressed as \6 — (5pred| + '^^pred < 0.004. Eor nuclei in the 
shape-transitional region this criterion is the last one to be satisfied among the four, because 
the potential energy surface (PES) is very flat versus 5. 

In order to accelerate the convergence when the wavefunction looks slowly converging to 
a state with deformation (5pred) we exert an external mass quadrupole potential during several 
tens of evolution steps so that the wavefunction quickly acquires the predicted deformation. 

^ In this paper, Spied is determined by an equation 5=0, where S is the time derivative of 5 and expressed 
as a quadratic function in 6 whose coefficients are determined from the least-square fitting to the last n steps. 
Changing n between 200 and f 50 produces a set of values {(5pred}- The average of them is adopted as Spred, while 
a half of the difference between the maximum and the minimum of them is used as the size of the error A5prcd- 
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Then, we switch off the external potential and continue the free imaginary-time evolution. With 
this acceleration method, the necessary number of evolution steps to achieve convergence can be 
decreased by a factor larger than two. The reason why this method works is that the quadrupole 
deformation is almost always the softest mode in atomic nuclei (under the symmetry). On 
the other hand, in generic multi-dimensional minimization problems, the principal difficulty is 



in finding the direction of the softest mode [ 31 1. 



In the top portion of Fig. || we show histories of the imaginary-time evolution of the defor- 
mation parameter 5 with (solid curve) and without (dot curve) the acceleration method for the 
prolate solution of ^IgErgg. The sharp fall of the solid curve is due to the external quadrupole 
potential exerted between the time steps 213 and 291. One can see that the convergence can be 
achieved within much shorter time than without the acceleration method. 

In the bottom-left and bottom-right portions of the figure, the time evolution of the max- 
imum spreading in the single-particle levels, Eq. (^), and the total energy (measured from the 
convergent value) are displayed, respectively, for the same process as in the top portion. One can 
see that the convergence of these quantities are also advanced by using the external potential. 



Figure |T] 



All the solutions of our extensive HF+BCS calculation have been found axially symmetric. 
The stability of these solutions against triaxial deformation can be known from the time evolution 
of the (very small) triaxiality parameter when the external potential is not effective, because it 
should grow exponentially if the axial solution is a saddle point rather than a minimum. It has 
turned out that none of the axial solutions exhibit clearly such an exponential growth of the 
triaxiality. It is possible, however, that triaxial minima exist which are separated from the axial 
path by a potential barrier. 

Bonche et al. found triaxial ground-state solutions for ®^Zr and ^'^Zr using the same Skyrme 
interaction and the same computer code as we used. However, these triaxial minima look so 
shallow that they can easily be moved to axially symmetric shapes by, e.g., making the pairing 
correlation slightly stronger. 

The treatment of the Coulomb interaction is in line with the appendix C of Ref. |3^: The 
exchange part of the interaction is approximated by a local potential in terms of the Slater 
approximation. The direct part is implemented in terms of a Coulomb potential obtained by 
solving the Poisson equation using a three-point approximation to the second derivative. The 
boundary conditions are determined by, first, dividing the nucleus into two fragments in one of 
the x-y, y-z, and z-x planes (we select the plane which gives the largest distance between the 
centers of mass of the two fragments) and, second, expanding the Coulomb potential originating 
in each fragment in terms of multipole moments of protons through order three. 

Throughout this paper, we treat a nucleon as a point particle, neglecting its extension. The 
effect of the spurious center of mass motion is corrected for by reducing the bare nucleon mass 
by a factor 1 — A~^. The explicit form of the hamiltonian density is found in Ref. 0]. 



3 The nuclear masses 

In the manner described in section ^, we have solved the HF-I-BCS equation for even-even nuclei 
with 2 < Z < 114. For each isotope chain, the calculation extends from outside the proton drip 
line by several neutrons to beyond the experimental frontier in the neutron-rich side by a few 
neutrons. For Z > 100, the calculation extends to < 2Z — 42. Spatially localized solutions 
have been obtained for 1029 nuclei. Although some of these solutions have negative two-proton 
separation energies, their densities are localized owing to the Coulomb barrier. The second local 
minima have been found for 758 nuclei. 
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We determine the ground-state solution of each nucleus by, first, searching the spherical, a 
prolate, and an oblate solutions and, second, comparing the energies of thus obtained solutions. 
Our strategy to search for these three solutions for each nucleus is as follows. The spherical 
solution is obtained by constraining the mass quadrupole moments to be zero. The prolate 
(oblate) solution is searched in two steps. First, we exert an external potential proportional 
to Qz on the initial wavefunction until its quadrupole deformation parameter satisfies 6 > Q.l 
(< —0.1). Second, we switch off the external potential, let the wavefunction evolve by itself (or 
with the acceleration method described in section 2.4), and see if it converges to a deformed local 



minimum. If the nuclear shape becomes very close to the sphericity in the course of evolution, 
i.e. 5 < 0.02 (> —0.02), we conclude that the normal-deformation prolate (oblate) solution does 
not exist in this nucleus. 

For some nuclei with 28 < Z,N < 50, the FRDM[j^ predicts very large deformations 6 ~ 0.4. 
In order not to miss such large-deformation solutions, we have done additional searches for all the 
nuclei in this region, in which we continue to exert the quadrupole potential until 6 becomes > 0.4 
(< —0.3) before starting the free evolution for the prolate (oblate) solution. These additional 
searches indeed produced large-deformation solutions. However, none of them are the ground 
states unlike in the results of the FRDM. 

In shape-transitional nuclei, the PES often has more than three normal-deformation minima. 
However, they are usually very shallow and it is doubtful that each of them corresponds to a 
distinct eigenstate notwithstanding the quantum fluctuation in shape. Therefore, we do not 
manage to find out all of these shallow minima. 

3.1 Correction for the finite mesh size 



As we have discussed in section 2.1, it is necessary to correct the total energy for the inaccuracy 
due to the finite mesh size because its relative error has to be by far smaller than that of other 
quantities for the sake of comparison with experimental data. 

To evaluate the size of the error, one needs the HF-I-BCS solution for vanishingly small mesh 
size. It can be obtained without difficulty concerning spherical solutions, because even personal 
computers can execute spherical HF-I-BCS codes with very small radial-grid spacing. Therefore, 
we have chosen to compare the total energies between the solutions of a spherical HF-I-BCS 
cod^ (the SKHAFO taken from Ref. [^9|) and the Cartesian-mesh code EV8 with constraints 
of vanishing quadrupole moments. In the following, we designate the total energy obtained with 
the spherical code as Eq while that from the Cartesian-mesh code simply as E. Our aim is to 
construct a formula for the energy correction so that the corrected energy Ec = E'-l- (correction) 
has a much smaller r.m.s. deviation from Eq than E has. 

In the least-square fitting to determine the parameters of the formula, we have used 1005 
nuclei (among the 1029 nuclei) whose paring gaps (both for proton and neutron) coincide within 
0.1 MeV between the results of the two codes. With the simplest fitting function of -Ec = ^ + 
ciA, the r.m.s. difference between Ec and Eq can be decreased from 6.7 MeV to 0.35 MeV. By 
adding temrs up to the second order in and Z, 

Ec = E + ciA + C2{N -Z) + c-iA^ + c^A{N - Z) ^ c^{N - Zf, (7) 

the r.m.s. value of E^ — Eq is decreased to 114 keV with ci = —40.2, C2 = —20.6, C3 = 0.033, 
C4 = —0.081, and C5 = —0.080 (keV). Since this size of error is much smaller than the typical 

^ The spherical code has been modified such that the following points are the same as in the Cartesian code: 
the values of physical constants, the treatment of the pairing correlations, and the modification of the proton 
potential outside the Coulomb barrier. The volume of the spherical cavity was set to the same size as that of the 
rectangular box used in the Cartesian code. The radial-grid spacing was set at 0.1 fm. 
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precision of the mass formulae in the market place (~ 0.5 MeV)||^, we have decided to 

adopt the above formula^. 

We have tested the accuracy of the correction formula (|^) for deformed solutions by com- 
paring Ec with an energy extrapolated to a ^ 0. In the extrapolation, first, we calculate 
the total energy for seven values of a ranging from 1 fm to 0.56 fm (six values ranging from 
1 fm to 0.6 fm for ^^''Pu), while keeping the box size constant. Second, we fit a polynomial 
E{a) = E'ext + bio^ + 62a® to the seven or six sets of values (a, E) using Ee^t, bi, and 62 as the 
fitting parameters. This form of the fitting function has been chosen on the following grounds: 
At a ~ 1 fm, the error is dominated by a term of order a®, which originates in the seven-point 
approximation to the first derivatives. At a ~ 0.5 fm, the contribution from lower-order terms 
becomes comparable to that of the term. These terms seem to arise principally from the 
error in the Coulomb energy, whose leading order term is a^. The ambiguity in the extrapolated 
energy E'ext is roughly estimated to be ~ 0.2 MeV. 

In Table |l], the difference between the extrapolated energy to a — > and the energy calcu- 
lated with a=l fm is shown for oblate, spherical (obtained with constraints of vanishing mass 
quadrupole moments), and prolate solutions of five nuclei. The fifth column shows the energy 
correction given by formula (^). For these 15 solutions, the mean and the r.m.s. values of the 
difference between the correction formula and the extrapolations are 0.0 MeV and 0.21 MeV, 
respectively, which are smaller than or of the same size as the accuracy of the extrapolation. 
This agreement shows the applicability of the correction formula (^) to deformed solutions as 
well as to spherical ones. It also confirms that the spherical code and the Cartesian code with 
the spherical constraint are indeed equivalent to each other. 

Table 1 



3.2 Comparison of the masses 

Fig. ^ presents the corrected total energies Ec of 1029 even-even nuclei calculated with the 
HF+BCS method with the Skyrme SIII force^. For graphical reason, the smooth (i.e., macro- 
scopic) part -Emacro has been subtracted, which is defined by fitting functions of the Bethe- 
Weizsacker type, 

^macro = ayA + ag^'/^ + ai{N - ZfA-^ + acZ^A-^/\ (8) 

to Ec separately for yl < 50 and yl > 50, varying oy, 03, ai, and ac as free fitting parameters. 
The solid (open) circles designate that the nucleus is more (less) bound than the macroscopic 
trend -Bmacro, while the diameter of each circle is proportional to \Ec — -Emacrol- 

We also show the two-proton (two-neutron) drip lines for the HF+BCS with SIII and for the 
FRDM||l|. The two-proton (two- neutron) drip line lies between two even-even nuclei whose two- 
proton (two-neutron) separation energies 52p (»S'2n) have different signs. As for the two-neutron 
drip line for the SIII force, we use the macroscopic energy (Eq. (^) fitted to the SIII results (ay 
=-14.702 MeV, as = 14.05 MeV, ai = 21.47 MeV, ac = 0.6554 MeV) because our calculation 
does not extend to the neutron drip line. 

^ Inclusion of higher order terms of A'' and Z to the correction formula does not substantially decrease the 
r.m.s. error. With terms of degrees from zero to three (10 terms in total), the r.m.s. error is 107 keV. Addition 
of up to sixth-order terms (28 terms) leads to an r.m.s. error of 93 keV. On the other hand, addition of only one 
more term C6-Etdj4^/^ decreases the r.m.s. error to 86 keV, where JStd is the space integral of BspAp + BapnApn 
+ BeppApp (see Ref. for the definitions of B5 and Be). Further addition of a term C7(An + A.p)A reduces the 
error to 77 keV, where An and Ap are the pairing gaps. 

^ The nuclear binding energy corresponds to — -Ec, the nuclear mass to Ec + ZrUp + Nrrin, and the atomic 
mass to Zrric — Bc{Z) added by the nuclear mass. The values of nip, nin, and me are taken from Ref. ] 6|. We use 
Bc{Z)=U.33 Z^-^^ eV. The definitions of energy-related quantities like Sip and 5211 are given in Ref. |p4[ . 
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One can see regions of enhanced stability around double-magic nuclei with (N, Z) = (50, 50), 
(82, 50), and (126, 82). Another double-magic nucleus (82, 82) is outside the two-proton drip line. 
The super-heavy double- magic nucleus (184, 114) does not look like a local minimum of nuclear 
mass in this result. 

The two-proton drip lines of the HF-I-BCS with SIII (solid line) and the FRDM (dash line) 
are overlapping in most places. The distance between them is AZ = 4 for = 40, AA = 4 for 
Z =42 and 78, and AZ < 2 and AA^ < 2 for the other isotope and isotone chains. 

The two-neutron drip lines of the two theoretical approaches are also close to each other. The 
difference looks of the same size as that between the FRDM and the TUYY mass formula [p^, 
both of which are models whose parameters were determined through extensive fittings to the 
mass data. This fact indicates the quantitative appropriateness of the macroscopic isospin 
dependence of the SIII force. Indeed, it is what we expected in choosing the SIII force for our 
first extensive HF-I-BCS calculation. 



Figure || 



The r.m.s. deviation of the calculated ground-state masses from the experimental ones of 
480 even-even nuclei (the best recommended values of Ref. excluding those estimated from 
systematic trends) is 2.2 MeV. Note that the inaccuracy of calculations due to the finite mesh 
size remaining after the correction is by far smaller than this deviation. 

The difference for each nucleus is shown in Fig. |3|. The solid (open) circles are put when 
the calculated masses are smaller (larger) than the experimental ones, while the diameter of 
the circles is proportional to the magnitude of the difference. One can see that the calculated 
masses tend to be overbinding for Z=8 and 20 isotopes, and A=50, 82, and 126 isotones. Unlike 
spherical nuclei including these semi-magic isotopes and isotones, deformed nuclei have positive 
errors, which are ~ 3 MeV rather independently of the size of deformation. 



Figure |3| 



In Table g we show the r.m.s. differences of the masses of even-even nuclei between theoretical 
models as well as between the experiments and the models. In the table, AW'93 represents 
the experimental atomic mass table by Audi and Wapstrap], TUYY the mass formula of 



Tachibana et al. |33], FRDM the finite-range droplet model[y|, ETFSI the extended Thomas- 
Fermi Strutinsky integral method with the SkSC4 forceg, EV8C the HF-^BCS results using the 
Skyrme SIII force with the correction (Eq. (0)), and macro the Bethe-Weizsacker type function 
fitted to AW'93 (ay =-15.280 MeV, as = 16.01 MeV, ai = 22.33 MeV, ac = 0.6896 MeV). 
In parentheses are the number of nuclei to calculate the difference. The r.m.s. deviation from 
AW'93 is 2.2 MeV for EV8C, which is 3-4 times as large as that of 0.52 MeV for TUYY, 0.68 
MeV for FRDM, and 0.74 MeV for ETFSI. It should be noticed, however, that the parameters of 
FRDM, TUYY, and ETFSI were fitted to all the available recent experimental mass data while 
the parameters of the SIII force were determined by fitting to the masses and charge radii of 
only seven spherical nuclei. In addition, the number of the fitting parameters is 275 in TUYY, 
19 in FRDM, and 8 in ETFSI, while it is only 6 in the SIII force. 



Table 2 



As a further investigation of the macroscopic properties of the mass, we examine the possi- 
bility to decrease the r.m.s. deviation by improving the macroscopic part of the mass models. 
Namely, for each combination of the nuclear mass models and the experimental data, we add 
the Bethe-Weizsacker type function (Eq. (^)) to one of them and determine the four coefficients 
to minimize the r.m.s. difference between them. The resulting r.m.s. differences are tabulated 
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in Table |3|. This simple correction method can decrease the r.m.s. error of EV8C from experi- 
ments by 27%. As for the other models, however, the improvements are marginal. On the other 
hand, the differences between the models are greatly decreased because they come predomi- 
nantly from nuclei near the neutron drip line. This fact suggests that substantial improvements 
of the macroscopic part of the nuclear mass models are possible only if new experimental mass 
data of neutron-rich nuclei are provided. 

Table | 



3.3 Pairing gaps 

Let us compare the experimental and the calculated pairing gaps. The former are calculated by 



applying the 5-point formula|35] to the experimental masses^] of nuclei except those at major- 
shell closures (of protons for the proton gap, of neutrons for the neutron gap), while the latter 
are obtained directly from the BCS part of the HF-I-BCS scheme. 

When the gaps are plotted versus the mass number A, the experimental ones fall close to the 
empirical formula A = 12 A~^^'^ MeV. The calculated gaps also fall near the empirical curve for 
heavy nuclei but they are raised for 50 < ^ < 100. In lighter nuclei {A < 50), the calculated gaps 
are located below the curve, which is due to our setting the maximum pairing force strength to 
0.6 MeV. A likely origin of the overestimation for 50 < ^ < 100 nuclei is that the pairing active 
space above the Fermi level (see Eqs. (^)) is smaller than the wavelength of the shell oscillation 
{hujosc ~ 41j4~^/^ MeV) for small A; such a situation does not fit to the continuous spectrum 
approximation. In future calculations for A < 100, it is desirable to improve the method to 
determine the pairing force strengths. 



4 Deformation 



Since we impose the reflection symmetries in x-y, y-z, z-x-planes (the symmetry) on our 
HF-I-BCS solutions, as explained in section 2.1, our solutions have only multipole moments with 
even angular momentum / and even magnetic quantum number m. 



4.1 quadrupole moments 

In Fig. 1^, we compare the magnitudes of the intrinsic electric quadrupole moments between our 
HF-I-BCS solutions and the experimental values deduced from B(E2)|, i.e., the reduced electric 
quadrupole transition probability from the ground state to the first excited 2"*" state 0. This 
moment is defined for the HF-I-BCS solutions as 

Qo= E {2z'-x'-y'), (9) 

protons 

where z-axis is the symmetry axis (All of our solutions have practically axial shapes, as described 
in section ^). The deduction of IQol from the B(E2)| is based on the rigid rotor model [0, 

For nuclei with Qq > 3.5 b, the agreement with experiment is excellent. The even-even nuclei 
having the largest intrinsic quadrupole moment is ^9gCfi54 (indicated by letter A in the figure) in 
the experimental table|^, while it is foe^^fiss among the 1029 HF-I-BCS solutions (Qo=16.6 b). 

The largest discrepancy is found in ^gQThi32 (indicated by D), whose experimental Qq is 5.5 b 
while the HF solution is spherical. Two nuclei located at isolated points are ^^gPtgg (indicated 
by B) and ^88Rai34 (indicated by C). 

For nuclei with smaller Qq, however, many nuclei falls not in the diagonal line but in a 
horizontal line, which means that the HF solution has a spherical shape when the experimental 
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B(E2)| is not necessarily very small. This may be explained by attributing the enhanced B(E2)| 
not only to static deformations but also to the collective shape oscillation around the spherical 
equilibrium. It may also be related to the complicated landscapes of the potential energy curves 
of nuclei with ^4=50-100 (see section 4.3). 



Figure H 



4.2 Deformation parameters 

For certain purposes, the deformation parameters are more useful than the electric multipole 
moments, although the former are model-dependent while the latter are directly related to 
the experimental observables. In a widely-used method of the Strutinsky shell correction, the 
deformation parameters are built in the theory in order to specify the single-particle potential. 
On the other hand, for mean-field solutions, one has to define the deformation parameters from 
the density distributions of nucleons. 

In this paper, we define the deformation parameters as those of a sharp-surface uniform- 
density liquid drop which has the same mass moments as the HF-I-BCS solution has. The mass 
density of the liquid drop is expressed as 

p{r) = poe{R{r)-\r\), (10) 
R{r) = RQ(l + Y,ai^Yi^{r)). (11) 

l,m 

The necessary and sufficient conditions on aim to fulfill the reality of R{r) and the D2h symmetry 
are that I and m are even numbers and a/^ = ofm ~ o^i-m- We set aim = for / > 6 and 
determine the remaining seven parameters po; -^O; ^20) 0,22^ ^40; 0^42; and 044 such that the liquid 
drop has the same particle number, mean-square mass radius, and mass quadrupole {r'^Y2m) 
and hexadecapole {r^Y^rn) moments^ as the HF-I-BCS solution has. 

The resulting values of Rq and 020 coincide very well with (|)^''^rrms and 5, respectively, 
where r^^^ is the r.m.s. radius and 5 is the deformation parameter defined in section For 5361 
samples (mass, proton, and neutron moments of 1029 ground and 758 first-excited solutions), the 
maximum and the r.m.s. deviations of Rq from (|)^''^rrms are 0.3 fm and 0.06 fm, respectively. 
Those of 020 from a fitted function (i|7r)V2^ -0.4752 -^0.785^ is 0.05 and 0.007, respectively. 

The axial quadrupole deformation parameter 020 is shown in Fig. |5| for the ground-state 
solutions of the HF-I-BCS equation with the SIII force. The open (solid) circles designate 
prolate (oblate) nuclei, while the diameter of the circles is proportional to the magnitude of 
the deformation parameter. The two-proton and two-neutron drip lines from our HF-I-BCS 
calculations (same as in Fig. |^ are drawn for the sake of convenience. 

One can see that nuclei at major-shell closures are spherical except for some N=28 isotones, 
while the deformation develops between the major-shell closures. Such pattern of the develop- 
ment of 020 looks regular for A > 100, while it is not so regular for A < 100, i.e., the change of 
deformation along isotope or isotone chains is not always smooth. In addition, oblate ground 
states are embedded here and there in light-mass region while for heavier-mass nuclei they are 
found only exceptionally in regions near major-shell closures. 

The largest deviation from the experimental 020 deduced from B(E2)| occurs in ^^C. For 
this nucleus the experimental B(E2)| is very large (corresponding to |a2o|=0.59) and an oblate 

^ One might wonder that the radial dependence of so strongly emphasizes the contributions from the 
peripheral regions that the deformation is sensibly affected by the shape of the box. To see the size of such 
erroneous effects we have done calculations by changing the radial dependence of the hexadecapole moment from 
r"* to in the determination of the deformation parameters. The resulting values of aim are only marginally 
changed except in a few nuclei outside the proton drip line. Incidentally, the FORTRAN source code to determine 
the seven liquid-drop parameters is obtainable via computer network. See an explanation at the end of this paper. 
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intrinsic deformation with a triangular three-alpha-cluster configuration has been suggested. In- 



deed, calculations using the Nilsson model[|36[ and the Strutinsky method |371| give oblate ground 
states. On the other hand, the HF+BCS calculation with the SIII force gives a potential energy 
curve which has only one minimum at the spherical shape. Other widely-used Skyrme forces of 
the SkM* and the SGII also give the only minimum at sphericity, while an old Skyrme force SII 



gives an oblate minimum with 6 = — 0.27[38]. However, one cannot deduce the superiority of 
the SII force since the optimal shapes of light nuclei are apt to be changed when effects beyond 
mean-field approximations are taken into account. For example, the parity projection may be 
important 1^] because the triangular three-alpha-cluster configuration violates the symmetry. 

A comparison with the results of the FRDM|jl| indicates a systematic differences in Z ~ ~ 
40 region, where our solutions tend to predict smaller deformations than those of the FRDM. 



We discuss on this difference in section 4.3, Another systematic difference occurs in a long and 
narrow region close to the proton drip line with 94 < Z < 102, where the FRDM predicts oblate 
shapes while our calculations give small prolate shapes, on which we do not discuss in this paper. 



Figure 5 



Deformations with the magnetic quantum number m=2 (triaxial deformations) are almost 
vanishing for all the nuclei we calculated: The magnitudes of 022 and 042 are smaller than 10~^. 
The deformation with m=4 (044) is larger than those with m=2. Its typical size is of the order of 
10-3 and the sign is mostly negative. However, these small non-axial deformations cannot affect 
any observables in experimentally detectable ways. Incidentally, it is an interesting problem how 
044 develops as the angular momentum increases [p9| . 

The axial hexadecapole deformation parameter 040 is shown in Fig. ^ in the similar manner 
as in Fig. |5[ This parameter becomes sizable for Z > 50. The sign of 040 is positive in the 
first half of the major shells and negative in the second half. This behavior is in agreement 
with the results of the FRDM|i| as well as with a naive expectation from the density profile of 
pure-j single-particle wavef unctions. The largest value of |a4o| is about 0.1 both for positive and 
negative signs. 



Figure |6| 



One of the advantages of mean-field methods over shell-correction schemes is that the protons 
and the neutrons do not have to possess the same radius and deformation. Making use of 
this advantage, we have calculated the liquid-drop shape parameters separately for protons and 
neutrons for 1029 ground and 758 first-excited solutions. As for 020, the r.m.s. difference between 
protons and neutrons is 0.012, while the maximum difference occurs in the ground state of ^iBeio 
(^20° ~ 0.52, O20" = 0.36). The absolute value of the difference is smaller than 0.05 for 98.9 % 
of the solutions. For 049, the r.m.s. difference is 0.0063, while the maximum difference is found 
in ^46612 (040° = —0.03, a^Q^ = —0.13). Concerning Rq, the r.m.s. difference is 0.18 fm. The 
maximum difference excluding nuclei outside the drip lines occurs again in ^^Be {Rq^° = 3.0 fm, 
i?Q^" = 3.9 fm). To summarize, the shapes of proton and neutron density distributions are not 
remarkably different from each other in any of our solutions. As for the radius, we also study 
the difference at much lower densities, i.e., skins and halos, in section ^. 



4.3 28 < < 50 nuclei 

Our calculations tend to predict smaller deformations than the FRDM|l|] for nuclei in 28 < 
Z,N < 50. These discrepancies in deformation often occur in a manner that our model gives a 
spherical shape while the FRDM predicts 6 ~ 0.4. These are originated in shape coexistence, 
i.e., the potential energy curve has more than one minimum which are energetically competing 
within ~ 1 MeV. 
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As an example to illustrate such a subtle situation, we take up 4QZr4o, for which the FRDM 
predicts a large prolate deformation a2o=0.43 while our calculation gives a moderate oblate 
deformation. In Fig. ^ the solid curve represents the potential energy curve calculated with the 
SIII force. It has as many as three minima, i.e., an oblate one at 6 = —0.18, almost spherical 
one, and a prolate one at 6 = 0.41. Because the energies of these three minima are so close to 
each other (they are within 0.6 MeV), the order of the energies can be altered easily by changing 
the parameters of the interaction. With our chosen parameters, the oblate minimum has the 
lowest energy (solid curve). If one decreases the pairing gap by 25% (i.e., the average pairing 
gap A,- with which to determine the pairing force strength is changed from 12 A^^/^ MeV to 9 
MeV) the prolate minimum becomes the ground state (dot curve). Instead, by changing 
the Skyrme force parameter set to the SkM*, while using the standard pairing force strength, 
one can make the prolate minimum the ground state (dash curve). On the contrary, the SGII 
force deepens the spherical minimum (dot-dash curve). The shapes of nuclei in this mass region 
have been studied in many papers |4^, 

Figure |^ 



4.4 Oblate solutions 

In Fig. ^, the energy differences between the oblate and the prolate solutions are plotted versus 
the neutron number. We put circles for those nuclei which have both the prolate and the oblate 
local minima. The circles belonging to the same isotope chain are connected by a line to guide 
the eyes. 

The energy difference is small near major-shell closures but large in the middle of the major 
shells. Apart from this shell fluctuation, one can see some overall trends. For each of the three 
major shells of neutrons divided by A^=50, 82, 126, and 184, the largest energy difference is 4.70 
MeV (^i|Sm62, 5=-0.25, 0.37), 6.58 MeV (i^|Gdioo,<5=-0.22, 0.31), and 11.05 MeV (fgtNoi52, 
5=— 0.17, 0.26), respectively; as the nucleus becomes heavier, the energy difference increases 
while the size of deformation decreases. 

For nuclei with N < 50, the changes along isotope chains are not so regular as in heavier 
nuclei and the oblate solutions often have lower energies than prolate ones. For nuclei with 
N > 50, oblate ground states are very rare and found only in nuclei very close to shell magics. 
The dominance of prolate deformations for > 50 may be attributed to the change of the 
nature of the major shells from the harmonic-oscillator shell to the Mayer-Jensen shell. 

Figure | 



5 nucleon skins 

As discussed in section ||, the greatest advantage of the mesh representation for this paper is 
that it is fit to describe skins and halos. 

We choose to define the skin by two conditions: A point r is in the proton or neutron skin 
when 

, . 3 , , , , , 0.16 fm~^ 
Prir) > -^Ptot(r) and ptot(r-) > — , (12) 

where r specifies proton or neutron and "tot" stands for proton-l-neutron. A similar definition 
was given in Ref. [^4|. The value of 0.16 fm^'^ is chosen as a standard density of the nuclear 
matter. Since the factors | and in the above definitions could have been differently chosen 
in rather arbitrary fashions, the absolute sizes of the skin thickness do not have a very precise 
meaning. Nevertheless, these values are useful to judge how the nucleon skins develop as the 
nucleus approaches to the drip line in the nuclear chart. 
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Fig. ^ presents the nucleon skin thicknesses calculated according to the above definition^. 
We put an open (solid) circle if the proton (neutron) skin exists for each nucleus. The diameter 
of the circle is proportional to the skin thickness. Among even-even nuclei inside the proton 
drip line, only nine have the proton skin (the heaviest one is 26Fe2o)- Proton skins thicker than 
1 fm are found in ^§04 (1.3 fm), 4662 (1.2 fm), and flSig (1.1 fm), all of which are single-magic 
spherical nuclei. On the other hand, the neutron skin has non-zero thickness in 39% of the 1029 
nuclei which we have calculated. In middle-weight nuclei ^5QSn84 and ^IgSnge, between which the 
present experimental neutron-rich frontier passes|45], the neutron skin is as thick as 1.02 fm and 



1.19 fm, respectively. The growth of the skin along isotope and isotone chains is monotonous 
and have no irregularity irrespective of proton or neutron and both inside and outside the proton 
drip line. 



Figure |9| 



When the nucleus is deformed, the thickness of the skin depends on the direction. (In Fig. ^, 
we take the arithmetic average of the thicknesses in the x-, y- and z-axes.) The ratio of the 
difference between the maximum and the minimum among the three axes to the average over 
the three axes is 8% for the proton skins and 20% for the neutron skins on the average. The 
relation between the anisotropy of the skin thickness and the quadrupole deformation parameter 



5 is shown in Fig. |10[ For each ground or first-excited solution having nucleon skins, a symbol 
is put at a point whose abscissa is 5 and its ordinate the skin thickness in the symmetry axis 
subtracted by that in the equatorial plane. The solid line is the result of a least-square fitting, 
while two dash lines correspond to twice as large mean-square deviation as the minimum value 
and can be used to judge the quality of the fitting. For neutron skins, there is a tendency that 
the skin is thicker in the symmetry axis than in the equatorial plane for oblate deformations 
and vice versa for prolate deformations: The neutron skin tends to make density more spherical. 
For protons, it is an interesting question whether the Coulomb repulsion between protons can 
reverse the trend in such a way that the proton skin promotes deformation. However, one cannot 



see any clear tendency in the right-hand side of Fig. 10, mainly because the number of points 
are fewer (heavy nuclei do not have proton skins). 



Figure |^ 



The nucleon halo, being composed of only two nucleons, has much lower density than the 
nucleon skin. We define the halo radius as the largest r satisfying the following condition, 

, , 0.16 fm-3 
Ptotir) > , 13) 

where ptot{f) is the angle-averaged mass density. Let us also define the halo thickness as the halo 
radius subtracted by 1.2 A^/'^ fm. According to these definitions, the halo radius (thickness) in 
the salient case of ^jjLig is 11.2 fm (8.5 fm) while the neutron skin thickness defined by Eqs. dlT 



is only 1.9 fm (from Fig. 4 of Ref. and Fig. 4 of Ref. I^^l)- -^^S- 11> show the relation 
between the skin and the halo thicknesses. One can see that the halo grows much slower than 
the skin for thin-skin nuclei (i.e. near the /^-stability line) but by far faster for thick-skin nuclei 
(i.e. near the drip lines). The maximum value of the halo thickness is 6.9 fm of 16S34 (indicated 

® In the mesh representation, the density is given only at the mesh points. In order to calculate the skin 
thickness in the a;-axis more accurately than the mesh spacing, first, we interpolate the density at the points 
of intersections between the a;-axis and the mesh planes x = (i + i)a using polynomials of order two in 
and 2^. Second, interpolations to arbitrary points in the x-axis are done for the logarithm of the density using 
polynomials of order three in x. The resulting density in logarithmic plots versus x shows no artificial ripples due 
to the interpolations. 
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by letter C in the figure) among the nuclei inside the drip lines. This value is somewhat smaller 
than 8.5 fm of ^^Li. One should note, however, that such a large halo as in ^^Li can be formed 
only in a very restricted interval of the last-occupied nucleon level [^]. It is possible that more 
gigantic halos are produced for different force parameters. 

Figure ^ 



6 Summary 

In this paper we have presented the results of our extensive Hartree-Fock+BCS calculations with 
the Skyrme SIII force for 1029 even-even nuclei with 2 < Z < 114 from outside the proton drip 
line by several neutrons to beyond the experimental neutron-rich frontier by a few neutrons. 

The single-particle wavefunctions are expressed in a three-dimensional Cartesian-mesh rep- 
resentation, whose advantages for atomic nuclei have been discussed. We have explained the 
newly developed points of the method of calculation, including the determination of the pairing 
force strengths, the acceleration of the convergence to the HF+BCS solution using an external 
quadrupole potential, and the correction of the total energy for the finite mesh size. 

We have compared the calculated ground-state masses with experimental data and with 
predictions by other mass models in section ^ The error from experiment is negative for most 
of spherical nuclei, while for deformed nuclei it is ~ 3 MeV rather independently of the size of 
deformation. The r.m.s. value of the error for 480 even-even nuclei is 2.2 MeV, which is much 
larger than the precision (~0.5 MeV) of recent extensively-parameter-fitted mass models like 
the TUYY and the FRDM. However, the proton drip line in the nuclear chart is almost identical 
to that of the FRDM. Even the location of the neutron drip line is not very distant from that 
of the FRDM and the TUYY. 

Deformations are discussed in section ^. Because of the D2h symmetry imposed on the 
solutions, only even electric multipole moments do not vanish. We have compared the elec- 
tric intrinsic axial quadrupole moments of our solutions with those deduced from experimental 
B(E2)|. The agreements between them are good except those nuclei with small quadruple mo- 
ments. After defining the deformation parameters aim for the HF+BCS solutions in terms of 
multipole moments, we have plotted the resulting 020 and 040 in the {N, Z) plane. The magni- 
tudes of non- axial deformation parameters 022, 042 , and 044 have turned out very small (|a22|5 
1 042 1 < 10~^, 1 044 1 ~ 10^^). The difference of shapes between protons and neutrons have been 
found small, too. Detailed discussions have been presented for ^^C and ^'^Zr. We have also 
shown the excitation energy between the oblate and the prolate solutions and pointed out a 
clear difference between below and above the A*" = 50 shell magic. 

Nucleon skins are discussed in section ^. We have shown that the skin grows monotonously 
and regularly as nucleons are added to the nucleus. On the other hand, the halo grows very 
slowly except near the drip lines, where it changes the behavior completely and expands very 
rapidly. It has also been found that the neutron skin tends to make the density distribution 
more spherical. 

Among the future problems for the extensive mean-field calculations of the nuclear ground- 
states properties are the improvements of the method to determine the pairing force strengths 
and the extension of the calculation to the neutron drip line by incorporating the coupling to 
the continuum in the pairing channel. Our final goal is the comparison of various Skyrme forces 
in their ability to reproduce the nuclear ground-state properties globally in the nuclear chart, 
and possibly the improvements of the force parameters. 

All the results of the extensive HF+BCS calculations reported in this paper are available 
electronically as an anonymous ftp service on the internet: ntl.c.u-tokyo. ac.jp. See a file 
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read.me in the home directory for instructions. The available quantities are the binding ener- 
gies, the Fermi levels, the pairing gaps, proton and neutron moments (r^, ^^12,01 ''^^2,2 '"^^,0! 
r'^Yi^2-, ?"^^,4), the skin thicknesses, the halo radii, and the proton/neutron/mass deformation 
parameters (020, 022, 040; 042, 044, Po) of 1029 ground states and 758 first-excited local 
minima. Single-particle spectra and density distributions of protons and neutrons are also ob- 
tainable for each solution. Additional postscript figures and some FORTRAN source codes to 
analyze the results are also provided. 
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TABLES 



Table 1. Estimated errors of the total energy due to the finite mesh size of a=l fm for oblate, 
spherical, and prolate solutions of five nuclei in MeV. The sign of the errors are inverted. 
The last column shows the values given by the correction formula See text for 
explanations. 





oblate 


spherical 


prolate 


correction 


^^Se38 


2.9 


3.1 


2.8 


2.83 


""^ 40 2^60 


4.4 


4.7 


4.2 


4.30 


'i^Ndro 


4.9 


5.1 


4.7 


4.99 




7.1 


7.3 


6.7 


7.14 


94PU146 


10.1 


10.0 


9.9 


10.05 



Table 2. The r.m.s. differences of nuclear masses between experiment and theoretical models 
and between theoretical models in MeV. In parentheses are the number of even-even 
nuclei to take the mean value. The FRDM (the ETFSI) does not give masses for < 8 
OT Z < 8 (for A < 36), while AW'93 and EV8C do not extend to the neutron drip line. 
See text for explanations. 



AW'93 TUYY FRDM ETFSI EV8C 



TUYY 


0.52 


(480) 












FRDM 


0.68 


(462) 


4.31 


(1521) 








ETFSI 


0.74 


(430) 


4.27 


(1472) 


2.74 


(1742) 




EV8C 


2.22 


(480) 


2.59 


(977) 


2.50 


(958) 2.26 


(940) 


macro 


3.55 


(480) 


17.25 


(2228) 


16.07 


(2246) 8.29 


(1895) 5.71 (1029) 



Table 3. Same as in Table |, but the r.m.s. errors have been decreased by adding Bethe- 
Weizsacker type functions to the models. See text for explanations. 





AW'93 


TUYY 


FRDM 


ETFSI 


TUYY 


0.50 








FRDM 


0.68 


1.73 






ETFSI 


0.68 


1.53 


2.10 




EV8C 


1.62 


1.80 


1.57 


1.54 



18 



FIGURE CAPTIONS 



Fig. 1. Imaginary-time evolution of the deformation parameter 5 (top portion), the maximum 
energy spreading of the single-particle states Ae (bottom-left portion), and the total en- 
ergy relative to the convergent value AE (bottom-right portion) for ^^^Er. The initial 
single-particle wavefunctions are taken from those of the HF+BCS solution for ^^^Er. The 
abscissae of the three graphs represent the number of steps of the imaginary-time evolu- 
tion. The dot curves represent the history of a free evolution while the solid curves are 
obtained by switching on an external quadrupole potential between time steps 213 and 291 
to accelerate the convergence. The quantity 6 is calculated in every step, while Ae and AE 
are in every 25 steps. 

Fig. 2. Nuclear masses calculated with the HF+BCS with SIII corrected according to Eq. (|^. The 
smooth part of Bethe-Weizsacker form has been subtracted. The solid staircase-like lines 
designate two-nucleon drip lines from our calculations. The dash lines are those from the 
FRDM|1|. 

Fig. 3. Error of the nuclear masses calculated with the HF+BCS with SIII corrected according to 
Eq. d?!). The grid indicates the locations of the magic numbers. 

Fig. 4. Comparison between the experimental and the calculated intrinsic quadrupole moments 
in unit of barn. For each nucleus whose B(E2)| is given in Ref. [0] (289 even-even nuclei 
ranging over 4 < Z < 98), a dot is put at a point whose abscissa is equal to the experimental 
value and its ordinate to the value calculated with the HF+BCS with SIII. The diagonal line 
is drawn so that one can see easily the quality of the agreement. See text for explanations. 

Fig. 5. The quadrupole deformation parameter a2o of the ground state solutions of the HF+BCS 
with SIII for even-even nuclei. Prolate (oblate) nuclei are designated with open (solid) 
circles whose diameter is proportional to |a2o|- The staircase-like lines represent two-nucleon 
drip lines from our calculations. 

Fig. 6. Same as in Fig. ^, but for the hexadecapole deformation parameter a^Q. The grid indicates 
the locations of the magic numbers. 

Fig. 7. The potential energy curves for ^''Zr obtained by solving the HF+BCS equation with con- 
straint on the mass quadrupole moment Qz- The abscissa is the deformation parameter 6, 
while the ordinate is the energy without correction for the finite-mesh-size effect. The solid 
and the dot curves are calculated with the SIII force, the former with the standard-strength 
and the latter with a weaker paring correlation. The dash curve is calculated with the SkM* 
force, while the dot-dash curve with the SCII force (vertically shifted by 15 MeV). 

Fig. 8. The energy difference between the oblate and the prolate solutions of the HF+BCS with 
SIII. The abscissa is the neutron number N, while the ordinate is the energy of the oblate 
solution subtracted by that of the prolate solution. A circle is put only when both solutions 
exist for each nucleus. The circles belonging to the same isotope chain are connected by a 
line to guide the eyes. The atomic number Z can be known from the direction of the hand 
in each circle. 

Fig. 9. The thickness of the proton and the neutron skins of the ground states of the HF+BCS with 
SIII. The proton (neutron) skin thickness is proportional to the diameter of open (solid) 
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circles. The staircase-like lines represent two-nucleon drip lines from our calculations. 



Fig. 10. The dependence of the anisotropy of the proton and the neutron skins on the quadrupolc 
deformation parameter 5. The circle (plus) symbol is used for nuclei with A < 100 (> 100). 
See text for explanations. 

Fig. 11. The relation between the thicknesses of halos and skins. For each of the ground and first- 
excited solutions of even-even nuclei, a plus symbol is plotted when the skin is of protons or 
there is no skins, while an x symbol is put when the skin is of neutrons. The plus symbols 
indicated by letters A and B correspond to fgSg and x4Si6, respectively, which are outside 
the proton drip line. The x symbol indicated by letter C corresponds to 16S34. 
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